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Abstract 

"Si In this work we estimate rates of the linear transient growth of the 

perturbations of cellular flames governed by the Sivashinsky equation. 
The possibility and significance of such a growth was indicated ear- 
lier in both computational and analytical investigations. Numerical 
investigation of the norm of the resolvent of the linear operator asso- 
^ ■ ciated with the Sivashinsky equation linearized in a neighbourhood of 

the steady coalescent pole solution was undertaken. The results are 
q | presented in the form of the pseudospectra and the lower bound of 

■ possible transient amplification. This amplification is strong enough 

to make the round-off errors visible in the numerical simulations in 
the form of small cusps appearing on the flame surface randomly in 
""^5 ■ time. Performance of available numerical approaches was compared 

to each other and the results are checked versus directly calculated 
norms of the evolution operator. 
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1 Introduction 



Sivashinsky's equation 

governs evolution of the perturbation $>(x,t) of the plane flame front moving 
in the direction orthogonal to the x-axis with the laminar flame speed u&, 
see Fig. ^ Here space coordinates are measured in units of the flame front 

oo 

width 5 t h, time is in units of S t h/ub, and H[&] = tt -1 / (x — y)~ 1 &(y,t)dy is 

— oo 

the Hilbert transform. 
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Figure 1: Perturbation (dashed line) of the plane flame front (solid line) 
moving at a speed u^. 

The equation was obtained in [l j considering the flame front as a surface 
separating combustible mixture of density p u and burnt gases of density p&. 
Assumptions of the low expansion rate pb/p u ~ 1 and small flame surface 
gradient |V$| <C 1 were also used in order to justify the appearance of the 
nonlinearity in where the parameter 7 = 1 — pb/p u - 

A wide class of periodic solutions to (JTJ) was obtained in 2J by using the 
pole decomposition technique. Namely, it was shown that 

$(x, t) = 27TNL- 1 (7 - 4vriVL- 1 ) t 

N 

+ 2 In |cosh [2vr& n (t)/L] - cos{2tt[x - a n (t)]/L}\ (2) 

n=l 

is an L-periodic solution to (JTJ if 

sin[27r(a TO - a m )/L] 

(it £ ^ \cosh[27r(6 n - b m )/L] - cos[27r(a n - a m )/L] 



da n 



2tt 



2 



+ 



sin[27r(a n - a m )/L] 



(3) 



cosh[27r(6 n + b m )/L] - cos[27r(a n - a m )/L] 



db. 



'n 



27TL- 1 coth (27ib n /L) - (7/2)sign6. 



dt 



a 




-± x \ cosh[27r(6„ -b m )/L)- cos[27r(a n - a m )/L) 
sinh[27r(6 n + b m )/L] 1 



sinh[27r(6 n - b m )/L) 



+ 



cosh[27r(6 n + b m )/L] - cos[27r(a n - a m )/L] 



(4) 



Here N is an arbitrary positive integer and prime in the symbol of summation 
means m ^ n. Pairs of real numbers (a n , b n ), n = 1, . . . , N are called poles 
and, correspondingly, function (j2J) is called iV-pole 1 solution to (JTJ). It is also 
convenient to consider <&(x,t) = const as a 0-pole solution to (fT|). 

If all the poles in (J3J), (@J) are steady and a n — a G R for n — 1, . . . , N, 
then, (J2J) is called a steady coalescent iV-pole solution. Solutions of the latter 
type, denoted here as $at(x) and illustrated in Fig. have been found to 
be the strongest attractors of ((H) and the period L preferred by ((TJ) has 
appeared to coincide with the size of the whole computational domain which 
we therefore denote as [— L/2, L/2], see e.g. 0. 

It was shown, see for example [4J, that for a given period L the num- 
ber of poles in steady coalescent pole solution (j2J) may not exceed Nl = 
ceil(7L/87r + 1/2) — 1, where ceil(a;) is the smallest integer greater or equal 
to x. Direct numerical simulations have revealed, in turn, that for suffi- 
ciently small values of L < L c the preferred number of poles is equal Nl. 
This observation was explained in [S] by means of the eigenvalue analysis 
of ((TJ) linearized in a neighbourhood of the steady coalescent pole solutions. 
The analysis has indicated that for any L > the steady coalescent iV^-pole 
solution is the only steady coalescent iV-pole solution to (0) with all the 
eigenvalues located in the left half of the complex plane. Strictly speaking, 
[3] does not provide a solid proof that their set of eigenvalues is complete 
and in this paper we explain why the comparison with the direct numerical 
calculation of the spectra, used in [Sj, cannot justify the completeness, in 
particular for large enough L. 

Surprisingly, for larger computational domains L > L c , numerical solu- 
tions to ((TJ) do not stabilize to any steady coalescent iV-pole solution at all. 
Instead, being essentially nonsteady, they remain very closely to the steady 
coalescent iV^-pole solution, developing on the surface of the flame front 

1 Strictly speaking, N is the number of complex conjugated pairs of poles a n ± ib n . 
However, we follow the tradition and keep this natural definition, as only real solutions 
are of interest. 
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Figure 2: Steady coalescent iV-pole solutions to the Sivashinsky equation. 
Here 7 = 0.8 and L = 1007T (N^ = 10). Graphs have been shifted vertically 
in order to get <fr N (±L/2) = 0. 



small cusps randomly in time, see e.g. (OJ. With time these small cusps move 
towards the trough of the flame front profile and disappear in it as this is 
shown in Fig. EJ 

The high sensitivity of pole solutions to certain perturbations was sug- 
gested in [7j as an explanation of the cardinal change in the behaviour of 
numerical solutions to (JIJ which takes place for L = L c . The argument of 
[Zj was based on a particular asymptotic solution to an approximation of 
the Sivashinsky equation linearized in a neighbourhood of the steady coa- 
lescent iV-pole solution. In the following works, see e.g. jH], the approach 
has been developed further and a model equation with stochastic right hand 
side, explicitly representing the noise, has been proposed and investigated. 
Sensitivity of Sivashinsky equation to the noise has also been studied in jO], 
and an estimation of dependence between L c and the amplitude of noise in 
the form of the the round-off errors has been obtained in JUj in a series of 
direct numerical simulations. 

Similar insufficiency of the eigenvalue analysis to interpret time depen- 
dent behaviour of asymptotically stable systems is also known from prob- 
lems of classic hydrodynamics, such as Poiseuille and Hagen-Poiseuille flows 
|10j . The failure of the spectral analysis in these problems was linked to the 
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nonorthogonality of the eigenfunctions of the associated linearized operators 
and was explained by the estimation of possible transient growth of pertur- 
bations [H], ^2]. A convenient tool to estimate possible transient growth 
of solutions governed by nonnormal operators were developed during the 
last decade in the form of the pseudospectra 13 . Corresponding numerical 
techniques have been reviewed in |14j . 

In this work we estimate rates of the linear transient growth of the pertur- 
bations of the steady coalescent iVVpole solutions to the Sivashinsky equa- 
tion. In Section 2 we linearize the equation in a neighbourhood of the steady 
coalescent pole solution. In Section 3, results of direct computations of the 
pseudospectra of the linear operator are presented. Also, a comparison of 
performance of available numerical techniques is given. Estimation of the 
rates of growth in terms of Kreiss constants, norms of the Co-semigroup and 
condition numbers, are presented in Section 4. We conclude with a discussion 
and a summary of results in Section 5. 



2 Linearized Sivashinsky equation 

Substituting t) = $>jsr(x,t) +<f>(x,t) into (JTJ) and neglecting terms which 
are nonlinear in <ft(x,t), one obtains 

^ = A N <j), t > 0, (5) 
where operator A n is defind by the following integro-differential expression 

A N u = q N — + — + -p- — tt - ^ xeR, (6) 

ox Ox z 2 ox 

on sufficiently smooth L-periodic functions with the square integrable on 
[-£/2,.L/2]. Here 

_d$ jV _47rJ^ sin[27r(x-a)/L] 
N dx L ~[ cosh(2vr6 n /L) - cos[27r(x - a)/L] ' 1 ' 

L > and a G R are real parameters, and the set b n , n = 1, . . . , iV is the 
steady solution to (jlj). 

The adjoint integro-differential expression is 

A > = -— + dx^ + 2^ (8) 

Hence, .4^*4.^ 7^ A* n An for > 0, and operator An is nonnormal. However, 
for 0-pole solution the first term in the right hand side of (JHj) disappears and 
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operator An is normal. In this sense we can say that it is the nonlinearity of 
the Sivashinsky equation what makes its associated linearized operator An 
nonnormal. 

If (JHJ), © is differentiated by x, then the resulting equation for tp(x) = 
d(f>/dx is dip /dt = A'n^P, where 

A N ip = — — + — T + -— ~ — , xeR. 

ox ox 2 ox 

The eigenvalue problem A'nv = Xv was studied in Obviously, eigenvalues 
of An and A'n are the same and the eigenfunctions of the latter one are just 
x-derivatives of the eigenfunctions of the operator An- 

In accordance with jH] , operator A'n has a zero eigenvalue associated with 
the x-shift invariance of The same is true for An as well. Moreover, zero 
is at least a double eigenvalue of An, because ((TJ) is also $- shift invariant. 
Here, x- and $-shift invariance means, that if <&(x,t) is a solution to (jlj). 
then, for any Ci, Ci G -R, function $(x + Ci, t) + Ci is its solution either. 

If only solutions with the period L are of interest, then they can be 

represented by the Fourier series </>(x, t) = E <}> k e i '** kx l L . Substituting 

fc=— 00 

these series into (J5J), © multiplying the result by e l2nnx ^ L , and integrating 
over the interval x G [— L/2, L/2], we obtain 

' H ' + ?*(*)+ £ m^v(fc-m)0 m (t), (9) 



dt \ L 2 L 1 7 L 



where ^(fc) = L _1 Xf^/ 2 ^N(x)e~ i27rkx ^ L dx and < 00. The integral can 
be written as a linear combination of the integrals of type Jq cos my (a — 
cosy^dy and the latter one was evaluated by using entry 2.5.16.33, p. 415 
of [T5] yielding 



N 

y N (k) = -z4vrL- 1 sign(fc)e- i2 ^ a/L £ e" 2 ^ |fc|/L , (10) 



n=l 



Introducing the^representation of (jSJ) in the Fourier space d<j)/dt = An4>, 
the Fourier image An of the operator .4.^ is defined by the (fc,m)-th entry 
of its double infinite (— 00 < k, m < 00) matrix as follows 

(Av)fc,m = (-4tt 2 L~ 2 A; 2 + tttL- 1 ^!) 5 kjTn 

N 

+ 87r 2 L- 2 msign(A; - m ) e -a<k-™)a/L ^ ^m^-HA | m | < ^ 

n=l 

(ID 
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where 5k )m is the Kronecker's symbol. 

It can be shown, that the value of the free parameter a does not affect 
neither spectral properties of An nor its 2-norms. Hence, we consider the 
case a = only. 



3 Pseudospectra of the linear operator 



3.1 Computational techniques 

In what follows we will work with matrix cut off at \k\, \m\ = K, i.e. 
all (k,m) entries of An with either \k\ or \m\ greater than K are neglected. 
Thus, instead of matrix An acting on double infinite vectors <ft, we consider 

the (2K + 1) x {2K + 1) matrix , whose entries coincide with those of 
A N for -K < k,n < K. 

In accordance with [Hj, in order to estimate possible nonmodal ampli- 



fication of solutions in (jSJ, we first calculate values of 



zl-A 



(K) 



A' 



as a function of the complex parameter z for large enough values of K. 

Level lines of this function form boundaries of the pseudospectra of An" \ 
see jEJ. It was found in numerical experiments that a good level of accu- 
racy of the most interesting part of the pseudospectra of A n is achieved if 
the cut off parameter K is about 2L/tt or greater. Thus, and in virtue of 



the Parseval identity ||(;zX— A 



(zT 



-A 



N 



, speaking about 



pseudospectra or other £ 2 - n orm based functionals of .Ajv we actually mean 

those calculated for An" ^ with large enough K. Here and in what follows 
1 is the unity operator or matrix of an appropriate size. 



Calculations of 



zl-A 



-{K) 



N 



can be carried out straightforwardly, 
however computational costs can be reduced if .4a/ ^ is transformed appro- 



priately. By construction, matrix .4^ acts on vectors u^' = (u-k 



Uq, Ui, 



Let us rearrange their components and consider w 



(K) 



[Bo, (s«) , (Sf») 



, where u± 



(u±i, . . . , u±k) t ■ The permuta- 



tion matrix V, corresponding to the proposed rearrangement w^ K ' = Vu^ K \ 

transforms An" ^ into VAn" ^V^ 1 , which acts on and has the following 
structure 



VA 



N 



V 1 




(12) 
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One zero eigenvalue of An , corresponding to the $-shift invariance of (|T|). 
can be seen from f)12j) explicitly. For other blocks of (fT2"|) we have: 

C m = -87r 2 L- 2 \m\ £ e -^b n \ m \/L^ m= i j 2,...,K, 

n=l 



V k , m = {-At: 2 L~ 2 k 2 + ^1-^)5^ k, m = 1, 2, . . . , K, 



A { ^ m = 8n 2 L~ 2 msign{k - m) £ e - 27rfe "l fc - m l/ L , fc, m = 1, 2, . . . , if, 



iV 



n=l 



.4^ = -87r 2 L- 2 msign(A; + m) £ e - 2 ^™l fc+m l/ L , jfe, m = 1, 2, . . . , if. 



A/" 



n=l 



Following the idea of [H], we apply the similarity transform 



/ 1 
X 

\ox 



(A") 

to PAv V~ l . Here X is the unity K xK matrix corresponding to the block 
structure of (|12jl . Unlike [Hj , the normalizing coefficient 1 / a/2 was chosen to 

preserve the 2-norm. The transformed matrix TV An {TV)- 1 is decoupled 
into two diagonal blocks 



TV An K) {TV)' 1 



/ C 

V + A^+A^ | , (13) 

V V + A {1) - 




and has the same 2-norm as An ■ The 2-norm of (|13|) is the maximum 
of 2-norms of its two blocks, each of which is of twice smaller size than 

TVAn^ \tV)~ 1 . In practice, the number of arithmetic operations re- 
quired to estimate the 2-norm of a matrix is of the order of the cube of 

(K) 

its size. Therefore, estimation of the 2-norm of An through blocks of 

TV An \TV)~ l is more efficient. From our experience, the 2-norms of the 
blocks are of the same order of magnitude, although the 2-norm of the upper 
block supersedes the lower one for most of practically important values of z. 

A straightforward and reliable way to calculate the 2-norm of the resol- 
vent of An* ^ (or of diagonal blocks of (fI3]) ) is through the singular value 
decomposition (SVD). Namely, the reciprocal to the smallest singular value 



8 



s of zl — An ^ is equal to 



zl-A 



-{K) 



N 



sec 



The direct Matlab 



implementation of SVD worked well in our case, though a few inverse iter- 
ations with continuation in z, suggested in [IB], appeared to be as accurate 
and, on average, about six times faster. 

An alternative algorithm is based on projection to the interesting sub- 
space through the Schur factorization followed by the Lanczos iterations. It 
was suggested in [T3j in the form of a Matlab script and is, on average, about 
two times faster than inverse iterations with continuation. Further, our tests 
have shown that its efficiency degrades much slower as the matrix size or 
required accuracy grows. Thus, Schur factorization with Lanczos iterations 
was the algorithm of our choice. It was intensively monitored by the direct 
SVD, however. 

A comparison of performance of the inverse iterations with continuation 
and of the Schur factorization with Lanczos iterations is given in Fig. El for 



calculations of 



zl-A 



(K) 



N 



with L = 40vr, 7 = 0.8 and K 



80. 



The criteria of stopping the iterations was 



» 



/4"» 



< 0.01, i.e. 



in) 



when the relative increment of the n-th approximation s 

(K) 

singular value Sq of zl — An is smaller than e = 0.01. 



to the smallest 



Graphs reveal 

areas with slower convergence of iterations. Unlike the number of required 
inverse iterations is usually less than the number of the Lanczos ones, the 
latter are much cheaper computationally, resulting in a significantly better 
overall performance. 



3.2 Structure of the pseudospectra 



As we are interested in stability of the steady coalescent pole solutions and 
possible rate of linear growth of their perturbations, the vicinity of the 
imaginary axis is of principle interest. Reflection symmetry of the func- 



tion 



[zl-A 



in regard to the real axis proves it sufficient, for our 



G {z : —5 < 



< 1, < 9f(z) < 1} 



,-1 



for L = 407T and 



purposes, to study it in the region z 
only. 

Figures andEt illustrate level lines of (zl — An) 
7 = 0.8. A rectangle plotted with a dashed line in Fig. marks the location 
of the area magnified in Fig. |3Jd. Asterisks in the figures show approxima- 
tions to the eigenvalues of the operator An- These parameters correspond 
to the appearance of microcusps in our direct numerical simulations with 
single accuracy. The picture suggests that the large area of high values of 



(zl - A 



N 1 



near the origin can be the reason of significant amplification 
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in 
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Figure 3: Relative increment of the smallest singular value in the last iter- 
ation (top) and number of iterations (bottom) carried out by methods (TOJ 
(left) and [H| (right). 



of the round-off errors, which in the case of single accuracy are of order 10~ 7 . 

Second critical case corresponding to the appearance of microcusps in 
calculations with double accuracy, see Fig. El is shown in Figs. Eh andEJi. 



{zT - A 



Nj 



near z = suggests a 



re round-off errors which are of order 



Again, a large region of huge values of 
possible match with the magnitude of t 
10~ 16 in this case. 

Data on the 2-norm of the resolvent of An for L = 2007T and 7 = 0.8 are 
given in Fig. 0^. The figure shows further widening of the area of large values 
near the real axis. Accordingly, calculated eigenvalues 



(zl-A 



Of 

spread further from the real axis and to the right from the imaginary axis. 



Also, they tend to form a cluster near the level line 



[zl-A 



N) 



10 



15 



10 




11 



cf I 

Also, Fig. E£ demonstrates that our calculations fail to produce reliable 

(K) 

results if the smallest singular value s of zX— An is less than about 1CT 15 . 
This is because of the effect of the round-off errors of the computer on the 
computational algorithm used to estimate sq. We see, however, that level 
lines of (zT — An) 1 2 corresponding to so > 10~ 15 are much less sensitive 
to these round-off errors than the eigenvalues. 

All the algorithms for estimation of pseudospectra mentioned in Section 
13. II are subject to the effect of the round-off errors and special arrangements 
are required in order to get reliable results for s < 10~ 15 . In particular, cal- 
culations with 128-bit arithmetic, implemented in some computer systems, 



'zl-A 



N 



corre- 



can be used. However, for our purposes knowledge of 
sponding to s > 10~ 15 is sufficient. 

The last example of the pseudospectra for L = 10007T shown in Fig. 

(K) 

EJ was calculated in a different way. The matrix An for K = 2000 was 
projected into its eigenspace spanning 1000 eigenvectors corresponding to the 
eigenvalues with the smallest absolute values. Then, the projected matrix was 



depicted in Fig. HJF. The 

oy asterisks, the neglected 



used to estimate the level lines of [zX — An) 
eigenvalues used for the projection are denoted 
ones by dots. The eigenvalue problem for the original matrix of size 4001 x 
4001 was solved by the Matlab implementation of Q-R-iterations. We also 
tried to apply Arnoldi iterations in accordance with [T7] , but could not make 
them convergent even for the projection subspaces of smaller dimensions and 
for smaller values of L. 

The direct calculation of the approximation to the spectrum of An, 

(K) 

namely eigenvalues of An presented in Figs. - EF with asterisks and 
dots, was undertaken by the Matlab implementation of Q-R-iterations. In 
the case of L = 90n six directly calculated eigenvalues are located to the 
right from the imaginary axis, see also jH] • The number of eigenvalues in the 
right half of the complex plane grows for larger L. However, the pseudospec- 
tra plotted in Figs. |Ui - Ef suggest that these unstable eigenvalues cannot 



be trusted. They appear in the vast area of large values of 



and, in accordance with e.g. 



[zl-A 



N 



,-1 



can be very sensitive to the perturbations 



as small as 10 
question. 



-16 



which is on the level of the machine zero in the case in 



Every eigenvalue A of A 



-(K) 



N 



eigenvectors of An 
by ka||£||2 at most 



can be associated with a condition number 
1 , where w and u are corresponding normalized left and right 

(K) m • '- r {K) 

(K) 



sec 



if A 



N 



. Then, eigenvalues of An will be perturbed 
is perturbed by matrix £ with small enough 
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||£||2- Figure El illustrates these condition numbers for L = 407r and 907T, 
making a very good match to the magnitude of perturbations of eigenvalues 
given in Fig. HJi. Note, the rightmost eigenvalues are worst conditioned. 



o L=40% 
L=90tz 




(K) 

Figure 5: Condition numbers of the eigenvalues A of An for 7 = 0.8 



(K) 

We would like to stress that because of the severe nonnormality of A n 
some of its directly calculated eigenvalues may have nothing in common with 
what they should be in absence of the round-off errors. A particular numeri- 
cal method can even worsen the estimation indeed. However, no one method 
can reduce the perturbation associated with the approximation of entries of 

An^ ^ by the finite-digit arithmetic of the computer, cf [Sj. The only way to 
increase the accuracy of the direct eigenvalue computations for L > 90n is 
to use a more accurate computer arithmetic with machine zero smaller than 

the reciprocal of max {^a}, where A(An^ "*) is the spectrum of .Aa/ \ 
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4 Estimation of the transient amplification 



4.1 Kreiss constants 



A robust lower bound on 



-Mn 



c 2 



can be obtained from the Laplace trans- 



form of e tAN , which under certain conditions (see JH]) can be written as 



J e- zt e tAN dt = (zl - A N ) 



-i 



Considering norms of both sides of this relation and carrying out straight- 
forward estimations of the integral: 



[zX-A 



N) 



\C 2 



J e~ zt e tAN dt 
o 



< sup He*- 4 ^!!/; 

t>0 



we arrive at sup \\e 
z with positive real part yielding 



c 2 o 
c 2 > 3l(z)\\(z2 — An)~ 1 \\c 2 - The latter is valid for all 



sup || e 



> 



sup 

5?(z)>0 



(z)\\(zl - A N )- l \\c 2 



Ani 



(14) 



where JC^ N is called the Kreiss constant, see also P^j . In other words, 
if JCa n is the Kreiss constant of the operator *4jv, then there is a pertur- 
bation (f)*(x,t) governed by (0) and a time instance £* such that the initial 
value of 4>*(x,t) is amplified at least JCa n times in terms of its £ 2 norm, i.e. 

||0*OM*)IU 2 > ^A N \\M x ^)\\c 2 - 

Our studies of pseudospectra represented, in particular, in Figs. Eh - EF 
indicate that the supremum in is reached on the real axis. Figure El 



shows dependence of the function z 



Z-zA 



(K) 



N 



on z for $s(z) 



0. 



We depicted results obtained by three different techniques and they are in 
good agreement with each other except for very small z. The discrepancy 
for z < 10~ 9 is because of the round-off errors as explained in the previous 

(K) 

should be of order 10~ 15 
1CT 9 . Indeed, this value of s 



section. The smallest singular value so of 1—zAn 

-l" 



to result in z 



1- zA 



-(JO 



N 



10 6 for z 



is too small to be accurately calculated on a computer with machine zero of 
order 10~ 16 . It is quite reliable to conclude in this case that fC^ N > 1.6 x 10 6 
for L = 407r and 7 = 0.8. 

Because of the effect of the round-off errors on the computation of so, 
similar estimations of the Kreiss constant for L > 80it on a computer with 



14 





\ 10 



2000 4000 
t 



versus z on the real axis 



Figure 6: An approximation to z [zl — An) 
for L = 40tt and = 160 (left). Dependence of £2 norms of Co- semigroups 

(K) (K) 

generated by An (solid lines) and En (dashed lines) on t (right). Here 
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0.8. 



the machine zero of order 10~ 16 are not accurate yielding a saturated value 
of order 10 13 . Instead, we have calculated more values of the Kreiss constant 
for a set of smaller L. Results are presented in Tableland Fig. El 



Table 1: Estimated Kreiss constants of An l 



L/n 
^A N 


10 

7.0 x 10° 


20 

3.3 x 10 2 


30 

2.1 x 10 4 


40 

1.6 x 10 6 


50 

1.3 x 10 8 


L/n 
^A N 


60 

1.2 x 10 10 


70 

4.0 x 10 11 


80 

5.0 x 10 12 


90 

1.5 x 10 13 


100 

1.5 x 10 13 



4.2 Norms of the Co-semigroup 



Good supplementary proof of essential nonmodal amplification can be pro- 
vided by direct estimation of the £ 2 norm of the Co-semigroup generated by 
An- Similar to the previous estimations, we have calculated the 2-norm of 

the Co-semigroup generated by An ■ Matlab's implementation of a scaling 
and squaring algorithm with a Pade approximation has been used in order to 

avoid calculation of the Jordan decomposition of An^ ■ The results revealed 
a good convergence for K > 2L/tt. 

As we have mentioned in Sectional operator An has a nontrivial null- 



space M(An)- Because of this 



-JAn 



does not decay for t — > 00 and, 
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moreover, it grows slowly because of the round-off errors. In order to re- 
move the effect of the null-space on the asymptotics of decay, and also, to 
demonstrate that the amplification observed in numerical experiments was 
not caused by that double zero eigenvalue, associated with the translational 

modes, we have projected An* into its eigenspace M (An^ ^ orthogonal 
to m(An {K) \. The 2-norms of the Co-semigroups generated by the result- 



ing operator, denoted here as Bn , are depicted in Fig. IHb alongside with 

(K) 

the similar data for the original operator An 

(K) 

Construction of Bn for larger values of L is complicated by difficulties 
with the accurate identification of the eigenfunctions corresponding to the 
zero eigenvalues. The latter ones appear to be perturbed and are as distant 

from z = as a few other eigenvalues. Note, that projection of An* ^ into 
M [An i affects the Co-semigroup not only asymptotically for t — ► oo, 



but for t — > as well. 

Data in Fig. |Hb matches our estimations of the lower bound of the pos- 
sible amplification of perturbations to the Sivashinsky equation and their 
extrapolations for larger values of L. Also, they show that presence of the 
nontrivial null space, corresponding to the shift invariance of the equation 
is not responsible for high sensitivity of the steady coalescent pole solutions 
to the noise. The latter conclusion is reinforced by the comparison of the 

(K) (K) 

pseudospectra of An and Bn ■ On scales of Figs. andEJo they are 
simply indistinguishable and can only be seen in a very close proximity of 
the origin, as shown in Fig. 

One may see that the only effect of the projection is a small shift of the 
pseudospectra to the left, resulting in the reduction of the Kreiss constant 
of about 30 times. It is still well above 10 4 , however, perfectly matching the 
corresponding curve in Fig. Eb- 



4.3 Condition numbers 

A traditional estimation of the Co-semigroup generated by is given by 



exp { t inf Mz)] } < 



c 2 



< K 2 (A N )explt sup [£(z)]>, (15) 



zeA(^jv) 



where A(^4at) is the spectrum of An, see [19J. If An is a finite-dimensional 
operator, then ^{An) is the condition number ^{An) = cond 2 (V) = 
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xlO~ 




xlO~ 



-0.02 
*(z) 



Fi gure 7: Comparison of level lines of log 10 

-l" 




log 



10 



zl-B 



-(K) 



N 



zl-A 



(K) 



N 



(left) and 



(right) for L = 40vr, 7 = 0.8, and K = 160. 



II V || 2 1| V^ 1 1| 2 of the matrix V whose columns are formed by the eigenvec- 
tors of A n- For infinite-dimensional operators, meaning of k 2 (An) is not so 
straightforward and, what is even more disappointing, it is often infinitely 
large, see [T3]. However, we try to estimate ^(An), because if successful it 



would give an estimation of the upper bound of 

for sup = as follows: 

zeA(A N ) 



£2 



following from (fTK|) 



C-2 



< K 2 (A 



Ni 



(16) 



Figure |H1 depicts graphs of k 2 ( A 



(K) 



cond 2 (V^' ) ) versus L for dif- 



r{K) 



ferent cut off parameters K. Here columns of matrix Vjy are eigenvectors 



-{K) 



-(JO 



of matrix An ■ The difference between k 2 \^An j for different K may 

look small on the graph. However, the graph is in the log 10 scale and the 
discrepancy is on the level of an order of magnitude. Hence, convergence is 

not obvious and we do not pose obtained k 2 \An I as an estimation of the 



upper bound of 



in (fTBJl . 



The graph of cond2 B 



(K) 



versus L is also illustrated in Fig [HI Unlike 



k 2 (An), which estimates the upper bound of amplification of solutions of the 
initial-value problem for the number cond 2 ( £>a/ j gives an estimation 
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10 20 30 40 50 60 70 80 90 



Figure 8: Dependence of condition numbers of V N and Bjy (marked with 
asterisks) on L for 7 = 0.8. 

of possible amplification of perturbations of the right hand side / in the 
solution Uf of the linear equation £>a/ \ = f. Note, because of x- and 
$-shift invariance of ([TJ, condition number of An itself is infinite. 

5 Comparison of the estimations 

It was established in numerical experiments, see e.g. [I], that for small 
enough computational domains of size L < L c numerical solutions of 
stabilize to the steady coalescent iV^-pole solutions of ((TJ). This observation 
is in the explicit agreement with the eigenvalue analysis of the linearized 
problem carried out in [Sj. 

For larger L > L c , numerical solutions do not stabilize to any steady 
solution at all. Instead, being essentially nonsteady, they remain very closely 
to the steady coalescent A^-pole solution, developing on the surface of the 
flame front small cusps randomly in time. With time these small cusps move 
towards the trough of the flame front profile and disappear in it as can be 
seen in Fig. El 

Numerous numerical experiments did not reveal any significant depen- 
dence of the critical length L c on parameters of the computational algorithm. 
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They have shown, however, that L c is effectively affected by the round-off 
errors [H] . Thus, if / is the order of the amplitude of perturbations associated 
with the round-off errors, then L c = L c (f). Two values of L c (f) obtained 
in our calculations with 32- and 64-bit arithmetic are shown in Fig. EH Am- 
plitude / of the perturbations was of the order of machine zeros, i.e. 10 -7 
and 10 -16 correspondingly. Note, that in calculations with 32-bit arithmetic 
round-off errors dominated discretization errors [2*Uj . 




ylVn 



Figure 9: Dependence of the variety of measures of the critical strength f c 
of perturbations on the flame size L. 

It is convenient to invert the relation L c = L c (f) and write it in the 
form f c = f c (L), where f c is a critical noise strength for given size L of the 
flame. Reciprocal of the Kreiss constant K,a n , obtained in Section l4~T| can 
be considered as the lower bound of this critical strength f c of perturbations 
for any particular value of L. Here, the strength of the perturbation means 
its 2-norm. Corresponding graph is plotted in Fig. El It is in a very good 
agreement with the results of our direct numerical simulations. The graph of 

^2 [A-n" ^ versus L is also given in Fig. El for the illustrative purposes. We 

remind, that there was no evidence of convergence of k 2 yA^ ^ to k 2 (An) 

in our calculations and interpretation of the graph as the upper bound (fTB*|) 
of f c is not justified. 
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An analytical attempt to estimate the value of f c was made in [7j where 
the following modification of ©, © has been considered: 

dcj) x d(f) d 2 (f) idH{4>] 

+ tj o + -—k — , x e R. (17) 



dt R N dx dx 2 2 dx 

Here Rn = (d 2 $N / 'dx 2 ) 1 is calculated in the crest of the steady coalescent 
A-pole solution, see Fig. |2J In [7j a particular asymptotic solution to (JT7J) 
has been investigated. As a result, the dependence between the critical value 
of curvature radius Rn in the crest of the flame profile and the spectral 
density p/ of the most dangerous harmonics of <f)(x, 0) has been obtained. A 
functional link between L and Rn l of the steady coalescent iV^-pole solution 
to the Sivashinsky equation can be easily established yielding 

p /iC = 4-Ve^^/ 8 . (18) 

Here C\ and C2 are coefficients of the least squares fitting of Rn l = Rn l {L) 
with a straight line C\L + C2. 

When comparing our results with estimation ()18|) . the following should 
be taken into account. First, relation ()18j) has been obtained for the spectral 
density of the most dangerous harmonics of the perturbation (p(x, 0) rather 
than for its amplitude /. Second, assumptions made to obtain fT%|) are better 
justified for large L. The last but not least factor is that ()18)) is based on 
a particular solution and is likely to produce an overestimated value of p/ jC 
rather than the optimal one. In view of these peculiarities, the agreement 
between ffTHj) . obtained in [7j, and our estimations is striking. 

In contrast, the estimation of f c obtained in [H] is obviously out of the 
harmony. That estimation was based on studies of the dynamics of poles 
governed by ©, ©• Namely, the amplitude of perturbations to the solutions 
of was linked to the 6-coordinate of poles in the (a, 6)-plane. Then, 
analysis of the dynamics of these noise generated poles yields the estimation 

f c = 2 n 7r 6 7 - 5 L- 6 . (19) 

There is no doubt that the sensitivity of system ©, ® to noise is totally 
different of what we have for Analysis of the Jacobian of the right hand 
sides of system (J3J), (HJ) for the steady coalescent A^-pole solution reveals that 
this is a symmetric matrix and there is no linear nonmodal amplification of 
noise in (jHJ), (jlj) at all. For small L, when the nonmodal amplification is 
not essential, estimation (|T9"j) is in a good agreement with other data indeed. 
However, for larger L, the discrepancy between the results of [5] and of others, 
clearly seen in Fig. El can be interpreted as the measure of the importance of 
the linear nonmodal amplification of perturbations in the Sivashinsky equa- 
tion. 
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6 Conclusions 



In this paper we have undertaken the numerical analysis of norms of the resol- 
vent of the linear operator associated with the Sivashinsky equation linearized 
in a neighbourhood of the steady coalescent pole solutions. Performance of 
available numerical techniques was compared to each other and the results 
are checked versus directly calculated norms of the evolution operator. 

The studies demonstrated the robustness of the approach by resolving the 
problem of stability of certain types of cellular flames. They showed that the 
round-off errors are the only effect relevant to the appearance of the micro 
cusps in computations of large enough flames. These essentially nonlinear 
micro cusps are generated through the huge linear nonmodal transitional 
amplification of the round-off errors. Their final appearance and dynamics 
on the flame surface is governed by essentially nonlinear mechanisms intrinsic 
to the Sivashinsky equation. 

In order to retain its physical meaning for large flames, Sivashinsky equa- 
tion should be refined by accounting for the physical noise, e.g. in the way 
suggested in jH]. 
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